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1.  INTRODUCTION 


During  the  three  years  of  this  University  Research  Initiative,  7  papers  have  been  published  in 
quality  scientific  journals  with  high  standards  of  review  and  3  have  been  submitted  for  publication 
based  on  full  or  partial  support  of  AFOSR  89-0462,  including  the  work  of  four  Ph.D.  candidates 
(one  of  whom  received  his  Ph.D.  degree  within  the  last  2  years).  The  following  list  also  includes  a 
report  (partially  supported  by  this  URI)  on  a  software  package  for  solving  a  general  second  order 
parameter  involved  system  of  elliptic  partial  differential  equations  on  a  two  dimensional  domain 
with  Dirichlet  and/or  Neuman  boundary  conditions. 
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Hele  Shaw  Convection  with  Imposed  Shear 
Part  I;  Boundary  Layer  Formulation 


Ruby  Krishnamurti 
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Florida  State  University 
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Huijun  Yang 
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Chicago,  IL  60637  USA 


Abstract 

There  is  a  general  interest  and  need  to  understand  the  interaction  of  convection  with  imposed 
shear  flows.  In  most  previous  studies  the  shear  extends  throughout  the  convecting  layer  with  the 
result  that  the  convection  takes  the  form  of  rolls  with  their  axes  aligned  with  the  imposed  flow. 
However,  it  is  the  case  when  the  convective  flow  has  vorticity  not  aligned  with  the  imposed  flow, 
that  interesting  momentum  transfer  is  expected  between  convective  and  imposed  flows.  In  a  well- 
developed  shear  flow  with  the  shear  confined  to  the  boundaries,  non-aligned  convection  appears 
to  occur  in  the  interior  in  regions  of  negligible  shear,  (as  with  cumulus  convection  in  the  Earth’s 
planetary  boundary  layer). 

The  present  study  is  one  in  which  convection  is  forced  by  the  boundaries  of  a  Hele- Shaw  cell 
to  be  aligned  perpendicular  to  an  imposed  shear  flow.  The  imposed  shear  flow  may  be  a  Conette 
flow  extending  throughout  the  convecting  layer,  or  one  confined  to  a  boundary,  depending  upon 
the  geometry  of  the  Hele-  Shaw  cell.  In  Part  I  we  examine  the  case  in  which  the  imposed  shear 
has  a  boundary  layer  structure  and  study  its  interaction  with  the  convecting  interior.  In  Part  II 
we  will  present  the  results  of  small  but  finite  amplitude  convection  with  an  imposed  shear. 
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III.  FURTHER  COMMENTS 


1.  Convection: 

a.  Turbulent  convection  in  a  horizontal  layer  between  fixed  rigid  boundaries  sometimes  sponta¬ 
neously  produces  a  persistent  large-scale  shear  flow  superposed  on  the  irregular  pattern  of 
rising  warm  and  descending  cold  plumes  or  blobs.  This  shear  flow  implies  that  the  mean  heat 
transport  across  the  layer  is  accompanied  by  a  mean  flux  of  horizontal  momentum.  It  has  been 
suggested  by  Prof.  R.  Krishnamurti,  who  discovered  the  large-scale  flow  in  her  convection  ex¬ 
periments  around  1980,  that  analogous  phenomena  may  sometimes  occur  in  the  atmospheric 
boundary  layer  and  that  the  associated  momentum  flux  may  play  a  significant  role  in  the  large 
scale  dynamics.  However  other  things  beside  convection  are  involved  in  atmospheric  flows,  and 
it  seems  likely  that  a  somewhat  more  realistic  analog  would  be  convection  in  a  layer  one  of 
whose  boundaries  is  moving  horizontally  with  respect  to  the  other:  thus  the  fluid  motion  is 
driven  mechanically  as  well  as  thermally,  with  accompanying  flu.xes  of  heat  and  momentum.  In 
an  attempt  to  provide  some  theoretical  background  for  experimental  studies  of  this  interaction, 
Dr.  Howard  tried  to  extend  methods  previously  used  for  ordinary  convection  to  find  bounds  for 
heat  and  momentum  fluxes  in  terms  of  measures  of  the  strength  of  the  thermal  and  mechanical 
driving.  This  investigation  showed  that  there  is  indeed  a  limited  region  in  the  plane  of  the 
heat  and  momentum  fluxed  (whose  boundary  depends  in  an  explicitly  determined  manner  on 
the  temperature  and  velocity  differences  between  the  upper  and  lower  plates)  outside  of  which 
the  heat  and  momentum  fluxes  cannot  lie.  These  results  and  their  derivation  are  presented  in 
the  paper  by  L.  Howard  listed  in  section  I. 

b.  During  the  grant  period  R.  Krishnamurti  and  H.  Yang  have  completed  a  study  Shaw 

convection  with  imposed  shear.  The  problem  was  chosen  so  that  the  cellular  convection  is  forced 
by  the  boundaries  to  have  its  rotation  axes  perpendicular  to  the  walls  and  so  perpendicular  to 
the  imposed  shear  flow.  They  have  studied  the  boundary  layer  nature  of  the  imposed  shear  flow 
and  its  interaction  with  the  cellular  convection  of  the  interior  flow.  They  have  delineated  the 
conditions  for  steady  and  for  oscillatory  finite  amplitude  convection;  computed  heat  flux  and 
counter-gradient  momentum  flux.  During  this  period  also,  H.  Yang  successfully  defended  his 
dissertation  and  has  taken  a  post-doctoral  position  at  the  University  of  Chicago.  Department 
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of  Geosciences. 


2.  Stability: 

Work  by  Dr.  L.  Howard  directed  toward  finding  simple  methods  to  determine  the  number  of 
unstable  modes  in  hydrodynamic  stability  problems,  referred  to  in  a  particular  Ccise  in  the  first 
progress  report,  has  proceeded  somewhat  further  in  the  direction  of  generality  insofar  as  parallel 
flows  of  arbitrary  profile  are  concerned,  but  in  the  case  of  non-parallel  flows  no  significant  results 
can  be  reported. 

In  the  parallel  case,  the  numerical  method  of  integration  around  the  cut  to  determine  the 
number  of  unstable  modes,  which  was  described  for  the  case  of  uniform  shear  with  constant  Brunt 
frequency  and  Ekman  damping  in  the  first  report,  has  been  formulated  for  arbitrary  profiles,  and 
can  provide,  with  quite  simple  numerical  work,  the  number  of  unstable  modes.  This  is  useful  as 
a  guide  for  directed  numerical  search  for  unstable  modes,  enabling  one  to  decide  when  the  search 
has  found  all  of  them. 

A  second  aspect  of  this  is  that  there  is  an  old  result  relating  the  number  of  unstable  modes  of  a 
parallel  homogeneous  flow  having  one  of  a  special  class  of  velocity  profiles  to  the  number  of  negative 
eigenvalues  of  a  certain  related  Sturm-Liouville  problem;  but  the  way  that  result  was  obtained  is 
not  applicable  except  for  that  special  class  of  profiles.  In  the  past  year  Dr.  Howard  hais  found  a  way 
similarly  to  relate  the  number  of  unstable  modes  for  an  arbitrary  profile  to  a  kind  of  generalization 
of  a  Sturm-Liouville  problem.  This  generalization  seems  to  have  some  mathematical  interest  on 
its  own,  though  a  broad  overall  picture  such  as  we  have  for  ordinary  Sturm-Liouville  problems  has 
hot  yet  been  developed.  These  results  have  been  described  in  an  invited  talk  at  a  meeting  of  the 
Mathematical  Association,  but  have  not  yet  been  prepared  for  publication. 

3.  Rotating  Flow: 

The  quasigeostrophic  potential  vorticity  equation  is  an  approximation  to  the  full  rotating 
fluid  equations  that  is  appropriate  when  the  Ekman  and  Rossby  numbers  are  small  and  the  depth 
is  nearly  uniform.  As  mentioned  in  the  second  report.  Dr.  Howard  has  been  investigating  the 
approximation  that  is  appropriate  when  the  depth  is  no  longer  nearly  uniform.  The  results,  which 
are  obtained  in  a  way  that  is  conceptually  simple  though  technically  some  what  intricate,  are  quite 
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interesting.  In  some  ways  the  mathematical  structure  is  actually  simpler  that  for  the  nearly  uniform 
depth  case,  but  it  is  of  a  less  familiar  sort.  The  general  formulation  now  needs  to  be  applied  to  some 
specific  examples,  and  Dr.  Howard  is  currently  trying  to  work  out  the  details  for  cases  relevant  to 
some  of  the  experiments  on  the  annulus  with  a  “mountain”  which  have  been  done  at  GFDI  and  to 
which  Prof.  Pfeffer  and  Mr.  Ding  are  attempting  to  apply  the  quasigeostrophic  potential  vorticity 
equation. 

4.  Adaptive  Mesh  Refinement  Research: 

A  general  and  modular  adaptive  mesh- refinement  software  package  aimed  at  solving  general 
nonlinear  elliptic  partial  differential  equations  on  a  2-dimensional  domain  with  general  Dirichlet 
and  Newman  boundary  conditions  designed  by  Bank  (1990).  The  adaptive  procedure  of  mesh- 
refinement  was  implemented  using  a  finite-element  triangulation  using  a  picewise  polynomial,  using 
an  hierchicaJ  basis  multigrid  iteration  for  solving  the  resulting  systems  of  linear  equations.  We 
implemented  the  software  package  PLTMG  at  FSU  and  developed  a  full  graphics  systems  using  the 
IRIS  4000  color  graphics  to  display  the  results. 

2.  The  Turkel-Zwas  explicit  large-time  step  scheme  was  applied  to  a  hemispheric  barotropic 
model  with  constraint  restoration  of  integral  invariants,  in  spherical  coordinates  (Navon  and  Yu, 

1991) .  A  full  code  including  specialized  Fourier  filtering  for  polar  regions,  a  17-point  Shapiro  low- 
pass  filtering  combined  with  a  Robert  filtering  for  the  leap-frog  integration  were  applied  and  fully 
documented. 

A  novel  constraint  restoration  method  due  to  Navon  (1987)  was  implemented  as  well  in  order 
to  ascertain  that  the  quadratic  integraJ  invariants  of  the  shallow-water  equations  remain  almost 
unchanged  for  long-term  integration.  The  full  code  was  tested  with  various  initial  conditions  and 
was  run  along  with  graphic  routines  and  found  to  be  robust  and  performing. 

3.  Domain  decomposition  and  parallel  processing  of  a  finite-element  model  of  the  shallow- 
water  equations  (paper  submitted  to  Computer  Methods  in  Applied  Mechanics  and  Engineering, 

1992) .  A  non-overlapping  domain  decomposition  technique  wcis  applied  to  a  two-stage  Numerov- 
Galerkin  finite-element  model  of  the  shallow- water  equations  over  a  limited  area  domain.  A  Schur- 
complement  matrix  formulation  was  used  along  with  a  modified  interface  matrix  approach  in  order 
to  handle  the  coupling  between  the  various  subdomains.  A  preconditioned  conjugate  gradient- 


square  accelerated  iterative  technique  was  used  in  order  to  solve  the  resulting  systems  of  non- 
symmetric  linear  equations.  A  class  of  robust  and  efficient  boundary-probe  preconditioners  was 
employed  in  order  to  accelerate  the  convergence  on  the  interfaces  between  the  various  subdomains. 
We  then  parallelized  various  stages  of  the  finite-element  solution  of  the  nonlinear  hyperbolic  pole’s 
describing  the  shallow-water  equations  model  and  implemented  the  parallel  code  on  the  4  processors 
of  the  CRAY-YMP/FSU  supercomputer  using  multitasking  techniques  in  a  dedicated  environment. 

A  modified  interface  matrix  algorithm  approach  was  developed  based  on  modified  incomplete 
LU  factorization  (MILU)  preconditioner  to  accelerate  the  CGS  iterative  nonsymmetric  solver.  A 
speed-up  of  almost  a  factor  of  4  was  obtained  -  the  ratio  becoming  close  to  4  for  higher-mesh 
resolutions.  This  work  generalizes  the  domain  decomposition  approach  to  nonlinear  hyperbolic 
equations  and  presents  new  approaches  for  preconditioning  the  interfaces  of  the  subdomains. 
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IMPLEMENTATION  OF  PLTMG  VERSION  6.0 
ON  SCRI4D  VAX  WITH  IRIS  4000  GRAPHICSt 


By 


Jian  Yu 

Mathematics  Department 
Florida  State  University 
Tallahassee, FL  32310 


1.  Introduction 

PLTMG  is  a  software  package  for  solving  a  general  second  order  parameter  involved  system  of  el¬ 
liptic  partial  differential  equations  on  a  two  dimensional  domain  with  Dirichlet  and/or  Neuman  boundary 
conditions.  It  was  originally  designed  by  Bank  and  his  colIegues(1979b,  1982b,  1985b,  1988,  1990a)  as  a 
tentative  program  to  study  the  theoretical  and  practical  aspects  of  the  multigrid  iterative  method,  adaptive 
grid  refinement  and  error  estimation  procedures,  and  their  interation.  The  adaptive  procedure  embodied 
in  PLTMG  is  based  on  the  ideas  of  Babuska  (1986a,  1986b,  1986c),  Rheinboldt,  and  their  co-worker  (1986, 
1987).  The  resulting  systems  of  nonlinear  equations  are  solved  by  a  combination  of  the  approximate  Newton 
iteration  (See  Bank  and  Rose  1981,  1982)  and  the  hierarchical  basis  multigrid  iteration  (See  Bank, Dupont 
and  Yserentant,  1988,  Yserentant,  1985,  1986a,  1986b,  1986c).  The  mesh  refinement  algorithms  and  data 
structures  used  in  PLTMG  originally  followed  the  ideas  of  Bank,  Sherman  and  Weiser  (1983).  The  a  poste¬ 
riori  error  estimation  procedure  used  for  the  adaptive  mesh  refinement  was  developed  by  Bank  and  Weiser 
(1986,  1985,  1981).  The  pseudo-arclength  continuation  procedure  of  PLTMG  results  from  the  joint  work  of 
Bank  and  Chan  (1986).  The  sparse  Gaussian  elimination  procedure  and  global  stiff  matrix  storage  used  for 
the  resulting  linear  sparse  system  of  solutions  was  the  joint  work  of  Bank  and  Smith  (1987).  The  deflation 
technique  used  for  dealing  with  instability  caused  by  singular  matrix  was  developed  by  Chan  (1984). 

The  purpose  of  this  report  is  to  present  the  implementaion  of  the  software  package  PLTMG  Edition 
8.0  on  the  Florida  State  University  (FSU)  SCRI4D  Vax  and  to  present  some  examples  of  PLTMG  run  on  the 
SCRI4D  VskX  with  the  IRIS  4000  color  graphics  system.  The  main  algorithms  embodied  in  the  subroutine 
PLTMG  are  explained  in  detail  along  with  the  relevant  theoretical  background  and  major  results.  Data 
structures  used  in  PLTMG  are  also  illustrated. 
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Abstract:  We  present  non-overlapping  domain  decomposition  techniques  applied 
to  a  two-stage  Numerov-Galerkin  finite  element  model  of  the  shallow  water  equa¬ 
tions  over  a  limited-area  domain.  The  Schur  complement  matrix  formulation  is 
employed  and  a  modified  interface  matrix  approach  is  proposed  to  handle  the  cou¬ 
pling  between  subdomains.  The  resulting  non-symmetric  Schur  complement  matri¬ 
ces,  modified  interface  matrices  as  well  as  the  subdomain  coefficient  matrices  aje 
solved  using  Preconditioned  Conjugate  Gradient  Squared  (PCGS)  non-symmetric 
iterative  solvers.  Various  stages  of  the  finite  element  solution  are  parallelized  and 
the  code  is  implemented  on  a  foxir  processor  CRAY  Y-MP  supercomputer  applying 
multitasking  techniques  in  a  dedicated  environment. 


1.  Introduction  and  motivation 

Recently  there  has  been  an  increase  in  research  activities  in  the  area  of  parallel 
computing  due  to  the  advent  and  growth  of  various  parallel  processing  architectures. 
The  domain  decomposition  approach  achieves  the  highest  level  of  parallelism  in 
the  numerical  solution  of  partial  differential  equations.  Specifically,  the  problem 
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ABSTRACT 


In  this  study,  using  a  perturbation  method  we  first  obtained  the  steady  finite  amplitude 
solutions  and  their  stabilities  in  Hele  Shaw  convection  without  imposed  shear  flows.  We 
found  that  the  signs  of  and  are  opposite  to  signs  of  and  cr^~\  respectively, 
whereeis  and  are  zero.  This  result  shows  that  the  stabe  finite  amplitude  convection 
will  vanish  when  the  Rayleigh  number  is  larger  than  a  critical  value.  The  bifurcation 
diagram  accordingly  was  presented,  which  showed  that  there  was  another  branch  of  stable 
solution  when  the  finite  amplitude  increcises  to  a  critical  value. 

Using  a  modified  perturbation  method  with  a  two  parameter  expansion  and  a  strained 
time  coordinate  (PLK  method),  we  successfully  obtained  the  oscillatory  finite  amplitude 


solutions  in  Hele  Shaw  convection  with  imposed  shear  flows.  The  two  parameters  are  e 


and  Pe,  where  e  is  the  amplitude  of  the  stream  function  to  leading  order,  and  Pe  is  the 
Peclet  number,  which  is  the  ratio  of  the  thermal  diffusion  time  to  the  time  for  advection 
by  the  imposed  shear  flow.  The  stability  of  the  finite  amplitude  solution  with  and  without 
shear  was  tested.  We  found  that  both  are  stable  to  infinitesimal  disturbances.  It  has  been 
found  that  in  all  shear  flows  Riq,  wio,  -Rii)  -^21  and  0^22  were  zero, 

where  and  u>rnn  are  the  coefficients  of  e”*Pe”,  respectively  in  the  Rayleigh  number 

and  the  frequency  expansions.  However,  a^oi)  -^201  -^021  ‘^21  and  R22  were  not  zero.  Their 
analytical  expressions  in  terms  of  the  Prandtl  number  and  A,  which  is  the  aispect  ratio, 
have  been  obtained,  wqi  and  were  the  frequencies  due  to  the  linear  and  nonlinear 

interactions  between  the  imposed  shear  flow  and  convection.  i?02  and  J?20  are  due  to 
the  linear  interaction  between  the  imposed  shear  flow  and  convection,  and  the  nonlinear 
convection  interaction.  We  found  that  the  nonlinear  interaction  between  the  imposed 
shear  flow  and  convection  would  slow  down  the  pattern  propagation,  namely  reducing 
the  frequency.  In  the  Part  II,  we  will  discuss  the  stream  function,  temperature  pattern, 
generated  mean  flow  and  momentum  flux,  and  their  variations  in  Hele  Shaw'  convection 
with  imposed  shear  flows. 
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ASYMMETRIC  EKMAN  DISSIPATION, 
SLOPING  BOUNT) ARIES  AND 
LINEAR  BAROCLINIC  INSTABILITY 

HENG-YI  WENG  and  ALBERT  BARCILON* 

Geophysical  Fluid  Dyncniics  Instiiuie.  Florida  Stare  L'niversny. 
Tallahassee.  FL  S2J06.  L'S.A 

fRi'Ceiicd  /-  Dt'ccinher  :990:  in  mwt  form  19  Ann/  199!) 

*h;  enicts  of  svmrr.:;r:c  2nc  asymmetric  Hitman  cissioa:;on  on  baroclinic  insiabiiliy.  nnase  sceed  and 
»avi  iiructure  in  a  linear  Eaay-liice  channei  mode]  aitn  wiihoui  opposiieiy  siopmc  nouncaries  are 
namined.  In  the  aosence  of  siopina  boundaries,  symmetric  Ei.man  oissipaiion  has  a  siaoiiizina  tendency 
10 f  ail  iaves.  When  :ne  Ekman  layers  are  asymmetric,  the  viscous  asymmetry  destabilizes  waves  b>  iii 
ctienainc  the  unstable  waveband  toward  both  Iona  and  short  waves,  and  lii)  increasing  their  growth 
dies,  compared  with  viscous  symmetric  case.  When  the  asymmetric  dissipation  is  very  small,  the 
cesiabtlization  .may  result  in  a  frictional  instaoiilty  for  short  waves  which  are  staoic  in  innsutij  case.  The 
iiscous  asymmetry  makes  waves  dispersive  and  their  structure  asymmetric  about  mid-depth.  However, 
j:  sense  of  the  viscous  asymmetry  and  the  sign  of  shear  do  not  affect  the  instability,  but  do  modify  ihe 
anection  of  phase  propagation  and  the  shape  of  the  'vave  structure. 

With  asymmetric  Ekman  dissipation  and  sioping  bpur.aaries.  frictional  insiabiiiiy  is  a  mixture  of  three 
T.echanismst  I'ii  symmetric  dissipation  in  the  presence  of  sioping  boundariest  (lii  viscous  asymmetry  tn 
f.t  absence  of  sioping  boundaries:  and  liiii  tne  sense  of  mscous  asymmetry  m  the  presence  of  sloping 
boundaries,  which  is  sensitive  to  the  wavenumoer  and  the  sign  of  shear.  The  slope  .-enoers  the  pr.asc  speec 
“ote  dispersive  dnd  sensitive  to  the  sign  of  shear  and  the  sense  of  the  viscous  asymmetr-. .  For  westerly 
'.-.sar.  the  slope  resuces  iincreasesi  the  wave  amplitude  anc  maxes  the  wave  .more  oaroctintc  marotropici 
i.v  'ower  lupperi  ievei.  For  a  given  wavenumber,  easttrtv  snear  uvnamics  .may  be  ceaucec  from  vvestcriy 
..:ar  ayn  imics  by  a  proper  change  of  the  sense  of  the  viscous  asymmetry. 

KEb  WORDS;  Barociinic  waves.  Ekman  dissipation,  slooing  aoundancs. 


I-  introduction 


^uuse  of  its  linear  and  nonlinear  analytical  simplicity,  the  Eady  (1949)  model  of 
'’srociinic  instability  has  served  as  a  "workhorse  '  to  test  various  aspects  of  that 
instability.  Yet.  we  feel  that  the  vanous  possible  alternatives  offered  by  this  modei, 
as  asymmetric  Ekman  dissipation  and  sloping  boundaries  in  particular,  could 
^’^nd  a  thorough  re-examination  which  we  offer  beiow  to  lay  the  foundation  before 
’’■^'nnding  with  nonlinear  studies. 

,  ^'sing  a  two-iayer  /I-piane  modei,  Hoiopainen  (1961)  considered  a  single  Ekman 
the  lower  boundary  and  found  that  the  effect  of  friction  reduces  the  amplitude 
'  perturbation,  except  within  two  narrow  wavebands  adjacent  to  the  two  invtsctd 
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The  weakly  .honlinear  evolution  of  a  free  baroclinic  wave  in  ;ne  presence  of  slightly  supercritical, 
scrucally  sheared  zonal  How  and  a  forced  stationary  wave  fteid  that  consists  of  a  single  zonal  scale  and 
an  arbitrary  number  of  meridional  harmonics  is  e.xamined  ■vithin  :hc  conte.xi  of  :he  conventional  two- 
layer  model.  The  presence  of  the  ipianetary-scalei  stationary  wave  introduces  zonal  variations  in  ihe 
iupercriiicalitv  and  is  shown  to  alter  the  growth  rate  and  asymptotic  eauiiibrium  of  the  isynoptic-scalei 
baroclinic  wave  via  two  distinct  mechanisms:  The  first  ts  due  to  the  direct  interaction  of  the  stationary 
wave  with  the  shorter  synoptic  wave  iwave-wave  mecnanismi,  ana  ;he  second  ;s  due  to  the  interaction 
v'f  the  synoptic  wave  with  that  portion  of  the  mean  lleid  that  .s  corrected  by  the  zonally  rectified 
stationary  wave  fiu.xes  (wave-mean  mecnanismi.  These  mechanisms  can  oppose  or  augment  each  other 
depending  on  the  amplitude  and  spatial  structure  of  '.he  stationary  wave  field.  If  the  stationary  '.vave 
field  IS  confined  primarily  to  the  upper  iloweri  layer  and  consists  of  only  the  gravest  cross-'iream 
mode,  conditions  are  favorable  (unfavorable!  (dr  nonzero  e'juilibrium  of  the  free  wave. 

In  addition  :o  the  time  dependent  heat  fiu.x  generateo  by  baroclinic  growth  of  the  free  wave,  its 
mteraction  with  a  stationarv  '.vave  field  consisting  of  i  wo  or  more  meridional  harmonics  generates  time 
dependent  heal  llu.ves  that  vary  with  period  .)f  ihe  free  wave.  However,  if  '.he  stationary  wave  field 
contains  several  meridionai  harmonics  of  suiliciently  large  amplitude,  the  :‘ree  baroclinic  wave  is 
destroyed. 

KEY  WORDS:  Baroclinic  instability,  jtaltonary  waves,  vvave-wave  and  wave-.mean  How  interactions. 


1.  INTRODUCTION 

■A  characteristic  signature  of  geophysical  iluid  systems  is  the  presence  of  several 
mutually  interacting  scales  of  motion  (e.g.,  Cotucci  et  ai..  1981;  Li  et  ai..  1986; 
PfefTer  et  ai..  1990).  Here  we  focus  on  the  Unite-amplitude  interactions  between  two 
scales  of  motion  that  characterize  the  large-scale  circulation  of  the  atmosphere;  the 
stationary  planetary-scales  and  the  transient  synoptic-scales. 
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abstract 

Complc-t  principal  component  analysis  is  applied  to  data  Troni  ilucc  laboraiory  cxpcnmenis  of  flow  O'er 
iwo-wavc  sinusoidal  lioltont  topography  in  a  thermally  driven,  rotating  annulus  of  iluid.  The  e.'tpcnments  arc 
conducted  at  the  same  imposed  temperature  contra.si  1ST)  and  at  three  diiTereni  rotation  rates  (!2).  In  each 
case,  the  intensity  of  the  '.vavc  activity  is  ma.ximum  downstream  of  the  two  topographic  ridges.  The  analysis, 
however,  reveals  a  fundamental  difTcrcnce  in  the  behavior  of  the  waves  at  lower  rotation  rates  than  at  the  highest 
rotation  rate.  At  the  lower  U's.  the  baroclinic  waves  travel  over  the  topographic  ndges  with  diminished  intensity 
and  amplify  on  the  other  side  of  each  ndge,  with  the  result  that  the  flows  downstream  of  the  two  ndges  are 
coherent.  At  the  largest  !!.  at  which  the  Rossby  number.  Ro.  is  very  small  and  the  friction  parameter. 
r  =  E''VRo  (where  £  is  proportional  to  the  Ekman  number),  is  rather  large,  the  waves  downst.'eam  of  each 
ndge  are  decoupled  from  those  downstream  of  the  other  ndge.  such  that  there  is  no  coherence  between  them. 
It  IS  thought  that  this  behavior  might  be  related  to  the  small  Rossby  radius  of  deformation  and  large  eiTective 
Skman  layer  dissipation  associated  with  baroclinic  waves  at  large  rotation  rates. 


1.  Introduction 

Comple.x  principal  component  (CPC)  analysis  can 
be  a  powerful  tool  in  revealing  behavioral  character¬ 
istics  of  otherwise  complicated  llows  in  the  atmosphere, 
the  oceans  and  other  geophysical  systems.  Discussions 
of  the  methodology  and  applicatiotis  to  the  atmosphere 
have  been  published  by  Wallace  and  Dickinson  ( 1972), 
Wallace  (1972),  Pratt  and  Wallace  (1976),  Barnett 
( 1983),  Horel  ( 1984)  and  others.  More  recently,  this 


*  Geophysical  Ruid  Dynamics  Institute  Contribution  Number  230. 
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form  of  analysis  has  been  used  by  Bemadet  ei  al.  ( 19S9 ) 
to  study  Dows  over  bottom  topography  in  laborator>' 
e.spcriments.  The  latter  study,  wliich  showed  that  the 
traveling  waves  were  modulated  on  the  scale  of  the 
topography,  sparked  the  first  author’s  interest  in  ap¬ 
plying  CPC  analysis  to  the  e.spenments  discussed  by 
Li  el  al.  ( 1986)  and  PfelTer  ei  ai.  (  1989).  The  present 
paper  is  devoted  to  the  analysis  of  three  of  these  e.\- 
periments,  which  were  performed  in  a  thermally  driven, 
rotating  annulus  of  fluid  with  two- wave  sinusoidal  bot¬ 
tom  topography  at  a  single  imposed  temperature  con¬ 
trast  (AT)  and  different  rotation  rates  (G).  As  in  the 
atmosphere,  one  observes  in  such  e-xperiments  a  to¬ 
pographically  forced  “planetary  scale”  wave  and  a  train 
of  smaller  “synoptic  scale”  waves  arising  from  the 
baroclinic  instability  of  the  How.  Accordingly,  it  seems 
worthwhile  to  pursue  such  experiments  as  a  means  of 
gaining  insight  into  the  contnbution  of  topographic 
barriers  (in  distinction  to  other  inlluences  such  as  dif- 
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abstract 

This  paper  presents  evidence  of  the  sensitivity  of  a  general  circulation  model  ( GC.M )  to  the  time-differencing 
scheme  employed  wnen  the  pnysicai  parametencations  and  space  discretization  are  not  changed.  For  this  purpose, 
two  time-marching  senemes — the  ieaofrog  ana  tne  .Vlatsuno  schemes — are  analyzed  and  tested  on  tne  National 
Aeronautics  and  Space  .Administration-Goddard  Laboratory  for  .Atmospnenc  Studies  (  .NaSa-GL.aS)  founn- 
order  GCM  m  terms  of  the  stabiiity  and  benavior  of  2-month-averaged  helds.  Linear  analysis  suggests  that 
Rossby  waves  are  slightly  damped  and  slightly  accelerated  when  the  .Matsuno  scheme  is  used  and  that  these 
effects  are  scale  selective,  being  smallest  for  the  longest  waves.  It  also  suggests  that  such  waves  are  accelerated 
less  and  are  not  damped  when  tne  leapfrog  seneme  is  used.  .An  empincal  onhogonal  function  analysis  of  the 
meridional  component  of  velocity  at  -Jb’.N.  keeping  at  least  70Fe  of  the  vanance.  .'eveais  less  shortwave  activity 
in  the  numerical  solution  with  the  Matsuno  seneme  but  does  not  lend  suppon  to  the  conclusion  i.hat  the  waves 
are  accelerated  less  in  the  solution  with  the  leapfrog  seneme. 

The  two-dimensional  Eiiassen-Paim  f£-P)  fiu.s  divergence  and  the  eddy-induced  mean  meridional  arcuiation 
are  found  to  be  stronger  in  the  simulation  wuh  the  leapfrog  time-differencing  scheme  than  in  the  one  with  tne 
■Matsuno  seneme.  suggesting  that  the  transient-wave  activity  is  damped  by  the  Matsuno  scheme.  On  the  other 
hand,  the  uhree -dimensional  sutionary-wave  ac.iMty  flu.x  m  the  .Northern  Hemisphere  simulateo  wim  the  .Matsuno 
scheme  is  more  intense  than  that  simulated  with  the  leapfrog  scheme,  indicating  that  the  stationary  waves  are 
more  rcoust  in  the  integration  wn.h  the  .Matsuno  scheme. 

The  GC.M  precipitation  when  iniegraied  with  the  leapfrog  sche.iie  is  much  more  intense  over  the  tropical 
western  Pacific  and  the  nonheastem  Pacific  and  less  intense  over  the  western  Nonh  Atlantic  Ocean.  The  kineuc 
energy  of  waves  with  wavenumber  greater  than  9  simulated  by  the  Matsuno  scheme  is  consistently  smaller  than 
that  obuined  bv  the  leapfrog  scheme.  These  results  give  evidence  that  climate  simulations  are  sensiuve  not  only 
to  physical  parametenzations  of  subgnd-scale  processes  but  also  to  the  numerical  methodology  employed. 


1.  Introduction 

General  circulation  models  (GCMs)  are  the  most 
elaborate  of  a  hierarchy  of  mathematical  models  used 
in  the  study  of  climate.  general  review  of  atmospheric 
GCMs  presented  by  Simmons  and  Bengtsson  ( 1984) 
emphasizes  the  importance  of  physical  processes  in  de¬ 
termining  the  behavior  of  climatic  systems.  GGMs  have 
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been  used  for  seasonal  simulation  (ShukJaetai.  1981 ) 
as  well  as  for  studying  climate  variability  on  ume  scaies 
of  a  month  and  upward.  In  panicular.  studies  have 
been  conducted  to  determine  the  sensitivity  of  atmo¬ 
spheric  GCMs  to  changes  in  physical  mechamsms  such 
as  surface  albedo  ( Chamey  1975),  sea-ice  limits  in  the 
.Arctic  (Herman  and  Johnson  1978).  low-frequency 
variability  (Chamey  and  Shukla  1981 ),  and  inadequate 
orographic  effects  (Wallace  et  al.  1983). 

Reviews  of  the  numencal  techniques  used  in  nu- 
mencal  weather  prediction  models  and  GCMs  have 
been  presented  by  Mesinger  and  Arakawa  ( 1976)  and 
Kasahara  ( 1979).  Numencal  expenments  m  which  the 
earth’s  climate  is  simulated  using  different  numencal 
schemes  can  play  a  valuable  role  in  clarifying  the  nature 
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EXSHALL;  A  TURKEL-ZWAS  EXPLICIT  LARGE 
TIME-STEP  FORTRAN  PROGRAM  FOR  SOLVING 
THE  SHALLOW-WATER  EQUATIONS 
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.Abstract — A  FORTRAN  computer  program  is  presented  and  documented  appi>ing  the  Turkel-Zwas 
explicit  large  time-step  scheme  to  a  hemispheric  barotropic  model  with  constraint  restoration  of  intezrai 
invariants  of  the  shallow-water  equations.  We  then  proceed  to  detail  the  alzonthms  embodied  in  the  code 
EXSHALL  in  this  paper,  particularly  algonthms  related  to  the  eificiency  and  stability  of  T-Z  scheme  and 
the  quadratic  constraint  restoration  method  which  is  based  on  a  variational  approach.  In  particular  we 
provide  details  about  the  high-laiiiude  liliering.  Shapiro  filtering,  and  Robert  filtenng  algonthms  used  in 
the  code.  We  explain  in  detail  the  various  subroutines  in  the  E.XSHALL  code  with  emphasis  on  algorithms 
implemented  in  the  code  and  present  the  flowcharts  of  some  major  subroutines.  Finally,  we  provide  a 
visual  example  illustrating  a  4-day  run  using  real  initial  data,  along  with  a  samole  pnniout  and  graphic 
isoline  contours  of  the  height  field  and  velocity  fields. 

Key  Words.  Shallow-water  equations.  Spherical  coordinates.  Explicit  .hnite-di.'ferences.  Constraint  restor¬ 
ation,  Filtenng  techniques. 


INTRODUCTION 

Recently,  a  considerable  amount  of  work  has  been 
dedicated  and  aimed  at  efficient  integration  of  shal¬ 
low-water  equations  in  view  of  using  these  methods 
in  numencal  weather  prediction  models.  In  order 
to  achieve  computational  accuracy  and  efficiency, 
most  methods  are  concerned  with  the  different  time¬ 
scale  of  the  advection  and  the  gravity-inertia  terms 
in  the  shallow-water  equations  model  separately. 
Semi-implicit  schemes  (Robert,  1979;  Burridge,  1975) 
and  split-explicit  schemes  (Magazenkov.  Shvets, 
and  Shneyerov.  1971;  Gadd.  1978a.  1978b)  are 
examples  of  those  methods.  In  the  split-explicit 
schemes,  a  substantial  computational  economy  is 
achieved  when  compared  to  usual  explicit  time 
integration  schemes. 

Turkel  and  Zwas  (1979)  proposed  a  space-splitting 
rather  than  a  time-splitting  algonthm  for  the  explicit 
integration  of  the  shallow-water  equations.  Their 
method  is  based  on  the  fact  that  the  fast  gravity- 
inertia  waves  contain  only  a  small  fraction  of  the 
total  available  energy  and  therefore  these  waves  can 
be  calculated  with  a  lower  accuracy  than  the  slow 
Rossby  waves,  that  is  on  a  coarser  mesh.  .An  appli¬ 
cation  of  the  T-Z  space  split-explicit  integration 
schemes  with  real  initial  data  is  presented  and  dis- 
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cussed  by  Navon  and  de  Viiliers  (1987)  and  us 
properties  are  discussed  further  :n  Neta  anc  Navon 
(1989).  .A  linear  transfer  function  analysis  of  the 
shallow-water  equations  in  spherical  coordinates  for 
the  Turkel-Zwas  e.xplicit  large  time-step  scheme  was 
earned  out  by  Neta.  Navon.  and  Yu  (19901. 

The  purpose  of  this  paper  is  to  present  a  practical 
FORTR.AN  code.  EXSHALL.  which  implements 
the  T-Z  explicit  large  time-step  scheme  for  the  shal¬ 
low-water  equations  in  spherical  coordinates  along 
with  constraint  restoration  methods  for  enforcing  a 
posierion  conservation  of  the  integral  invanants  of 
the  shallow-water  eouations.  The  computer  program 
is  explained  in  detail  in  connection  with  the  vanous 
algonthms  implemented  in  the  code  EXSH.ALL. 
This  paper  can  be  used  as  a  user's  guide  to  the 
program  EXSH.ALL  both  in  providing  a  bnef 
descnption  of  the  theory  as  well  as  detailed  program¬ 
ming  implementation. 

We  present  here  the  T-Z  scheme  for  the  shallow- 
water  equations  in  spherical  coordinates  and  its 
related  algonthmic  background.  'Vanous  filters  used 
with  T-Z  scheme  and  which  impact  on  its  stability 
aiso  are  presented;  a  detailed  descnption  of  the 
vanous  subroutines  m  the  code  EXSH.ALL  is  pre- 
jcnted;  and  a  typical  example  of  a  4-day  run  with  the 
program  EXSHALL  is  presented  along  with  graphi¬ 
cal  output.  Finally,  in  the  .Appendix  the  commented 
and  documented  FORTRAN  source  listing  the  code 
of  the  program  EXSHALL  is  attached. 
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1-  Introduction 

Suppose  a  fluid  whose  motion  is  adequately  described  by  the  Boussinesq  equations  occupies 
a  horizontally  infinite  layer  of  thickness  d.  The  relevant  physical  properties  of  the  fluid  are 
described  by  the  thermal  expansion  coefficient  a,  the  kinematic  viscosity  u  and  the  thermal 
diffusivity  «;  in  addition  there  is  the  gravitational  acceleration  g.  directed  downward.  We 
suppose  the  top  of  the  layer  (r  =  d)  is  a  rigid  motionless  plate  held  at  temperature  rero. 
while  the  bottom  (r  =  0)  is  held  at  temperature  AT  and  is  also  rigid,  but  may  be  m.ovins 
in  a  fixed  direction  in  its  plane  with  speed  W' .  Dim.ensionless  independent  variables  are 
introduced  by  using  d  as  the  unit  of  length,  and  the  thermal  time  d- / k  as  unit  of  tim.e. 
AT  is  chosen  as  unit  of  temperature,  and  n/d  as  unit  of  velocity,  as  is  commonlv  done 
in  convection  problems.  With  a  suitably  defined  dimensionless  pressure  the  dimensionless 
Boussinesq  equations  are  then 

cr~‘(U(  -r  u  •  V  u)  -h  VP  -  RTk  =  u 

V  •  u  =  0 

T-u--r=  v'T 

k  is  the  vertical  unit  vector  and  the  parameters  are  the  Prandil  number  a  and  the  Rayleigh 
number  R: 


la  I 
1  16) 
i  Ic) 


<7  =  i//«  I  2a) 

R  =  agd^Td^ /{kv)  ,26) 

The  boundary  conditions  are  taken  in  the  form 

r  =  0.  u  =  0  on  0  =  1  ,;3a) 

T  =  1.  u  =  C'i-|-X'^j  on  0  =  0  (36) 


in  which  U  and  V  are  constants  with  L'*  -P  V-  =  \V-  =  i  IF'd/^:)-:  thus  the  Reynolds 
num.ber  '^ased  on  the  speed  of  the  bottom  plate  is  Re  =  W'd/i/  =  W/cr.  W'e  shall  consider 
flows  -  solutions  of  the  above  equations  -  for  which  ail  relevant  horizontal  averages  (denoted 
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The  sisady.  hydrostatic,  invisctc.  Boussinesq  flow  of  j  stabiy  stratified  fluid  over  a  beil-ihaoeo  ridee  ;s 
investigated  'iiihin  the  framewor.v  of  a  linear  model.  The  three  layer  mooei  atmosphere  iniroauced  is 
such  that  the  Brunt-Vaisala  frequency  is  constant  in  each  laver  but  the  interfaces  of  the  middie  layer 
are  allowed  to  vary  gently  in  the  cross-ridge  direction.  In  essence,  the  mode!  can  be  tuned  in  botn 
virttcal  and  norizoniai  directions.  These  cross-ridge  \ariations  can  produce  significant  differences  in 
both  the  cross-ridge  surface  'vmd  and  the  surface  drag  compared  to  the  response  ooiainec  ny  use  of  a 
bonaonially  uniform  renecting  layer.  These  changes  are  sensitive  to  both  the  -erticai  location  of  the 
•middle  layer  and  to  the  slope  of  ns  lower  inierface  at  the  ridge  crest.  Many  of  these  features  are 
eitpiained  by  means  of  a  conventional  layered-model  anaivsis. 

KHh'  WORDS:  .Mountain  wases.  stratified  fiow.  Boussinesq  fiuid. 

1.  INTRODUCTION 

Typically.  Jhe  static  stability,  expressed  by  the  Brunt-Vaisala  frequency,  undergoes 
considerable  variability  in  the  vertical  direction.  Figure  1  shows  that  the  static 
stability,  S~,  e,xhibits  variability  over  scales  of  1-2  .km  in  the  troposphere.  Simiiar 
features  are  retained  between  the  two  stations,  as  in  the  region  between  6  and 
Skm  in  panel  b.  In  other  cases,  however,  the  apparent  layering  of  the  static 
stability  may  be  localized  in  the  vicinity  of  each  station,  and/or  the  same  features 
may  ascend  or  descend  in  altitude  between  GCT  and  DEN.  This  latter  possibility, 
which  may  characterize  the  features  in  panel  a.  is  examined  in  the  present  studv. 

We  consider  the  dynamics  of  hydrostatic  mountain  waves  in  the  presence  of 
static  stability  variations  in  the  cross-ndge  direction.  Such  states  will  consist  of 
regions  bounded  by  .v-dependent  interfaces  across  which  the  Scorer  parameter 
changes  abruptly,  remaining  constant  on  either  sides  of  these  interfaces.  These 
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